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Abstract 

A  numerical  method  for  solving  the  full  non-linear  Boltzmann 
equation  is  presentedj  and  applied  to  the  problem  of  shock  structure 
in  a  gas  of  elastic  spheres.   The  success  of  the  method  hinges  on 
the  systematic  use  of  Gaussian  quadrature  and  Hermite  interpolation. 


Ill 


Introduction .  The  Boltzmann  equation  describes  the  evolution  of  the 
one-particle  distribution  function  f  =  f(x,u,t)j  where  x,  with 
components  (x-,,  x„j  x_),  is  the  position  vector,  u  with  components 
(u-,.  Up,  u_),  is  the  velocity  vector,  and  t  is  the  time.   In  the 
case  of  a  gas  of  elastic  spheres  it  has  the  form 


(1)    II  +  (u.V^)f+  ^(F-Zn^f  =  -|  /|v.el(f^f;-ff')du'da3 

where  m  is  the  mass  of   a  particle,  o    its  radius,   V   denotes  the 
gradient  operator  with  respect  to  the  x  variables,  V   denotes  the 
gradient  operator  with  respect  to  the  u  variables,  F  is  the  external 
force,  e  is  a  unit  vector  pointing  in  the  direction  of  the  solid 
angle  element  dwj  V  =  u'  -  u,  a  bar  under  a  symbol  denotes  a  vector 
quantity,  and 


f  =  f(x,u,t) 


f  =  f(x,u',t) 
f+  '  f(x,u'*',t) 


where 


u   =  u  +  (V*e)e_ 
u    =  u  -  (V' e) e 


u  ,  u    are  the  velocities  before  collision  of  those  spheres  which 
after  collision  have  the  velocities  u  and  u ' .   Analogous  expressions 
can  be  written  for  other  kinds  of  Interpartlcle  force.   For  an 
elementary  discussion  of  this  equation,  see  [1^];  for  a  thorough 


discussion  see  e.g.  [2]  and  [5].   The  right  hand  side  of 
equation  (1)  will  be  called  the  collision  integral. 

It  is  the  purpose  of  this  paper  to  present  a  numerical 
algorithm  for  solving  Equation  (1)  and  to  apply  it  to  the  study 
of  the  structure  of  a  shock  in  one  space  dimension.   Generalizations 
of  this  method  to  problems  involving  very  strong  shocks,  more 
space  dimensions,  and  other  molecular  models,  will  be  also  dis- 
cussed.  It  will  be  seen  that  the  solution  of  the  shock  problem 
provides  a  key  to  the  solution  of  the  other  problems;  the  main 
difficulty  has  been  overcome  in  the  program  discussed  in  this 
paper.   Furthermore,  the  numerical  solution  provides  insight  into 
some  approximate  procedures,  in  particular  Grad ' s  thirteen  moment 
approximation  [6]  and  Mott-Smith's  bimodal  approximation  [12]. 

Unlike  the  work  presented  here,  most  previous  numerical  treat- 
ments of  the  Boltzmann  equation  relied  on  a  Monte-Carlo  technique; 
some  of  these  treatments  are  ingenious  and  interesting,  but  none 
can  be  considered  accurate.   See  [1],  [8],  [9],  and  [13]. 
Reference  [9]  is  particularly  helpful. 

For  any  function  (f)(x,u),  let  ^ix)    denote  the  integral 

(j)(x)  =  /(})(x_,u)  f  (x  ,u)d]i  . 
Some  of  the  quantities  of  interest  in  the  solution  of  the  Boltzmann 


equation  are  the  following  moments  of  f :   the  density  p(x)  =  1,  the 

—  1  ~2  — 

mean  velocity  u,  the  pressure  p  =  ^pw  ,  where  w  =  u  -  u,  the 

teinperature  T  =  p/pR,  where  R  is  the  universal  gas  constant,  the 


1   2 

pressure  tensor  p.  .  =Pw.w.,  and  the  heat  flux  vector  S  =  ^-pw  w 

Other  quantities  of  interest  are  the  Boltzmann  H  function 


H  =  /  f  log  f  du   , 

and  in  the  shock  wave  problem,   various  geometric  parameters  which 
characterize  the  shock. 

In  the  case  of  a  gas  of  elastic  spheres,  the  mean  free  path  is 


^  =  l/(  /2  7T  p  a^)  . 


We  shall  now  specialize  equation  (1)  to  a  form  appropriate  to  the 
shock  problem.   Let  T„  be  a  reference  temperature,  and  u„  a 
reference  thermal  velocity,  u„  =  /2RT„ .   Let  t  be  a  collision  time, 
"^  =  /  u^ ,  and  fi,  a  reference  density.   Introduce  the  non-dimensional 
variables 


X*  =  xA  ,  u*  =  u/Uq,  t*  =  t/x  ,  f*  =  u^f/pQ  ; 

substitute  them  into  equation  (1)  and  drop  the  stars.   Furthermore, 

pick  units  in  which  a  reference  mean  free  path  l/(  ^/2  i\      pL  a  )  is 

1.   Assume  F  =  0,  (no  external  forces),  and  allow  f  to  depend  only 

on  one  space  variable  x,  =  x,  and  two  velocity  variables  u.,  =  u 

and  u  ,  where  u  is  in  the  direction  of  x  and  u   is  in  a  direction  ortho- 
r  r 

gonal  to  u.   These  assumptions  imply  that  the  flow  is  invariant 


under  rotation  around  the  x-axls.   Under  these  assumptions 
f  =  f(x,u,u^,t)  satisfies 

(2)    if  +  u  ^  =  /^d^  sin  \>   /^  dX  /^"du'  /+~  du^(f,f:-ff')-|v-e|/2/2  TT 
where 

e  =  (cos  \) ,    sin  \)   cos  X,  sin  ^    sin  X), 

u'  =  (u'.u  cos  X,  u  sin  X). 

u  =  (u,  u^,  0), 

V  =  u'  -  u, 

u+  =  u  +  (V.e)e  -  (u+,u+,u+), 

U.+  '  =  u  -  (V-e)e  -  (u+'  ,  u+ '  ,  u+  ), 


u+  =  v4+2+u+2  , 
r     23 


r       2      3 
f  =  f(x,u',u^). 


f  =  f(x,u,u^). 


f_l_  =  f(x,u+,u+), 

f;  =  f(x,u+',uj') 


This  is  the  form  of  the  equation  we  shall  use  below,  although  the 
method  of  solution  applies  to  the  general  equation  (1)  as  well. 


Principle  of  the  method.  There  are  several  major  difficulties 
in  the  solution  of  (1)  o.r^  (2).   The  function  f  depends  on  a 
relatively  large  number  of  independent  variables-three  plus  time 
in  the  case  of  one-dimensional  flow,  six  plus  time   in  the  general 
case-so  that  if  (2)  is  replaced  by  a  system  of  algebraic  equations, 
their  number  will  be  large.   The  presence  of  the  fourfold  non- 
linear integral  ensures  that  the  algebraic  equations  will  be  not 
only  numerous,  but  also  very  cumbersome.   Efficient  computation  is 
clearly  needed.   Another  difficulty  stems  from  the  nature  of  the 
collision  term,  more  specifically,  from  the  integration  over  the 
angular  variables.   Suppose  f  is  represented  by  a  discrete  set  of 
values  assumed  on  a  discrete  set  Z  of  points  in  phase  space.   The 
integration  over  u  ,  u   becomes  a  sum  over  the  values  assumed  by 
f  on  Z.   The  integration  with  respect  to  0,X  becomes  a  sum   over  a 
discrete  set  0  of  values  of  e,X.   For  any  reasonable  choice  of 
Z  and  0,  the  arguments  of  f^.,  f|  will  include  points  not  in  Z. 
Thus   interpolation,  both  accurate  and  stable,  between  values  of 
f  on  Z,  will  be  required.   (Monte-Carlo  methods  avoid  the  inter- 
polation problem  at  a  heavy  price  in  accuracy.) 

These  difficulties  can  be  resolved  as  follows :   once  one 
resigns  oneself  to  the  need  for  interpolation,  there  is  no  need 
to  identify  Z  with  the  nodes  of  a  regular  mesh.   One  can  then 
evaluate  f  at  points  (x^,  u^  ,  u   .)  where  the  u.,  u,^  •  are  at  one's 
disposal.   In  particular,  one  can  choose  u. ,  u   .  to  be  the 
roots  of  a  polynomial  Pjt(u),  where  P   is  the  N-th  degree  member 


a  sequence  of  polynomials  P„ ,  P  ,...,P  ,...,  orthogonal  with 
respect  to  a  weight  W(u).   With  this  choice  of  u  ,  u   .  one 
can  interpolate  between  the  values  of  f  at  these  points  using 
the  orthogonal  polynomials  P  (x).   Such  Interpolation  is  stable; 
see  details  below.   Furthermore,  with  this  choice  of  integration 
points,  integrals  can  be  evaluated  by  an  appropriate  variant 
of  Gaussian  quadrature.   In  the  case  of  the  Boltzmann  equation, 
it  is  natural  to  use  the  Hermite  polynomials  H  (u)  given  by 

H„(u)  =  (-l)"c„  e^^  -^   e-^^  ,   c„  =  (2%!)"^''^ 


"  "     dun       '    n 


—  1  /? 
which  are  orthonormal  with  respect  to  the  weight  W(u)  =  it     e~^ 

i.e. 

tt"-^^^  /  H  (u)H  (u)e"^  du  =  6     ; 
n'  '  m'  '  n,m  ' 

-u^/2 
6    the  Kronecker  delta.   The  set  {H  (u)e     }  is  complete  in 
n,m  n 

L  (-ooj+oo).   See  [10]. 

The  preceding  remarks  lead  to  the  following  step-by-step 
procedure  for  solving  (2):   Let  At  be  the  time  step.   Assume 
that  at  time  t  =  nAt  f  is  given  by  a  series 

1      P   'l   "^2 
(3)       f(x,u,u^,nAt)  =  TT  ^(u^)-^  y—  ^  a..(x,t)- 

•H.((u-v'^)/u'^)H.(u  /u")exp(-((u-v^)^  +  u^)/(u'')^) 
1      xxjrx  X      r  X 


2 


n  n 

where  v^  is  the  center  of  the  expansion  and  u   is  its  scale.   The 

A  X 

n      n 
subscript  X  In  u  and  v  indicates  that  both  parameters  are  allowed 

X         X 

to  vary  with  x,  and  it  Is  assumed  that  a^^  •  =  0  for  i  _>  L,  ,  j  >_  Lp. 

Appropriate  v  ,  u  ,  L^ ,  L„  will  be  determined  below. 
X    x    1    2 

It  is  adequate  to  evaluate  a. .(x)  at  the  points  x  =  kAx,  k 
integer.  Ax  a  spatial  increment,  and  we  can  write 


a"^   =  a.  .  (kAx,nAt) 
Ijk    ^J 


The  density  at  time  nAt  is 


p'^(x)  =  aQQ(x,nAt)  =  a^^Cx)  , 


the  temperature  is 


(4)  t"(x)  =  (u'^)2((3/2)a2Q  +  2-1^2^2q  +  2^/2^g2  ^/SRp " 


and  the  mean  velocity  is 


(5)  u-^  =  v^  +  2-l/2u^a5o/p^ 


Our  aim  is  to  obtain  f(x,u,u  ,(n+l)At)  as  a  series  of  the  form  (3), 
but  possibly   with  a  new  scale  u^    and  a  new  center  v     .   To 

X  X 

n+1 
achieve  this  aim  we  evaluate  the  values  f.    of  f(x,u,u  ,(n+l)At) 

Ij  k  r 

at  the  points  x,  =  kAx,  u.  =  v'^"''-^  +  u^"^-^  C-,  u   .  =  u"^"*"^?  ,.  ,   where 
k       '   1     X      X    ^1'   r,j     X   ^j  ' 

E,.,E,.    are  roots  of  Hj^(u)  =  0.   The  value  of  N  remains  to  be  chosen. 


The  algorithm  for  evaluating  f^-u^  will  be  described  below. 

Given  f(x,u,u  ,(n+l)At),  the  coefficients  a^ .  are  defined  by 


(6)   aj"!"l(x)  =  7r-l(u^+l)-2  //  f  (x,u,u^,  (n+1)  At  )H.  ( (u-v5;+l)/u^'^l) 


H.(u  /u'^"^^)dudu 
J   r   X       r 


^-l(^n+l)-2  ;;    f  H.((u-v'^+l)/u'^+l)H  (u  /uf  ^) 
X  1  X    j   -^   X 


exp((u-v'^"^^)/u"+l)2  exp(ur/u!;^+l)^  . 


r/"x 


■exp(-(u-v"+l)/u"+^)^  exp  -  (u/u.^'^^)^   dudu^ 

X  r   A         r 


An  obvious  change  of  variables  reduces  the  last  integral  to  the 
form 

2  -u2 
//  g(u,u^)e  ^  e   ^  dudu^ 

which   can  be   evaluated  by   Gauss-Hermite   quadrature    (see    [15]), 
i.e.    using  a   formula 

2   -u^  N        N 

(7)      //   g(u,u^)e-"   e      ^dudu     =   YZJZ  g(C,,5.)w  w      , 
^  ^        i  =  0    j=0  1      J       ^    •J 

where  5.,  C  are  the  roots  of  H  (u)  =  0  and  w. ,  w.  are  appropriate 
weights.   Because  of  the  choice  of  quadrature  points,  g  is  already 


known  at  the  appropriate  points  {^.,E,-).      Formula  (7)  Is  exact 
If  f(x,u,u  )  has  a  Hermlte  expansion  of  the  form  (3)  with 
L  =  N,  L2  =  N  (see  [11]).   Thus  as  L  ,  L  are  Increased,  N 
should  be  Increased. 

There  remains  only  the  task  of  deciding  how  to  evaluate 
f. .,  .   This  Is  done  In  the  present  paper  by  an  explicit  formula 
of  the  form 


(8)  fjJJ  =  A  q.^   *   At  Q(f,f) 


where  A  Is  a  linear  operator  such  that  f"^   -  Af'^  approximates 


9t      d 
We  use 


-  u  ^  ,  and  Q(f,f)  is  an  approximation  to  the  collision  Integral 


Af.  .    =  (l-u^)f.  .  ,  +  u  ^  f .  .  ,  ,  .  V 
ijjjk       Ax'  i,j,k     AX   i,j,k+s(u) 


where  s(u)  =  1  If  u  <  0  and  s(u)  =  -1  If  u  >  0.   This  A  has  an 
obvious  initultlve  appeal,  but  since  It  Is  only  of  first  order 
accuracy  a  more  accurate  approximation  may  be  In  order  In  future 
work.   For  stability  we  must  of  course   have 


(9)  ^  maxlu.l  <  1. 

1   -■- 


In  the  evaluation  of  the  collision  Integral,  the  crucial  fact  Is 

that  f  Is  given  as  a  continuous  function  of  u  and  u  ,  and  thus   no 

r 

further  Interpolation  problems  will  arise.   Using  the  conservation 
of  energy  and  momentum,  the  representation   (3)  of  f,  and  an 


obvious  charge  of  variables,  the  collision  term  at  (x,u.,u   .) 
can  be  reduced  to  the  form 

-2  -u'2 


i'"r,j'  11  ^   -1 


C(x,u.,u   .)  /  d^   f^    dX  7"''°°  du'  /"'"°°  du '  G(  6  ,X  ,u' ,u '  )e~^   e 


where  C  it?  a  constant  which  depends  on  x,  u.  ,  u   ..  This  integral 
can  be  approximated  by  the  mixed  Gauss  and  Gauss-Hermite  quadrature 
formula 


(10) 


N^   N2   N^   N^^ 


1  =  0  j  =  0  k=0  £  =  0     1   J   i^   r,v,   1  J  K  X, 


where  the  9.,  x.  are  the  roots  of  the  Legendre  polynomials  of  degree 

N,  ,  N„;  u,  ,  u„  „  are  the  roots  of  H„  (u)  =  0  and  H^,  (u)  =  0,  and 
1    <:;    K    r ,  X,  jg  j\j 

W. ,  W.,  w  ,w  are  the  quadrature  weights,  see  [15]. 

-'-J     K   X. 

Our  method  can  thus  be  summarized  as  follows:  At  the  begin- 
ning of  each  step,   the  solution  f  is  given  as  a  Hermite  series 
of  the  form  (3);  f  at  the  next  level  is  evaluated  at  appropriate 
points  using  a  difference  scheme  for  the  linear  terms  and  weighted 
Gaussian  quadrature  for  the  collision  term.   The  new  values  are 
then  synthesized  into  a  Hermite  series.   A  variant  of  this  method, 
using  Monte-Carlo  quadrature  rather  than  Gaussian  quadrature  for 
the  collision  term,  was  presented  in  [3]- 

The  fundamental  difference  between  our  approach  and  Hermite 
expansion  methods  such  as  Grad's  thirteen  moment  approximation  [6] 
lies  in  the  fact  that  the  number  of  Hermite  polynomials  used  is 
not  fixed  in  advance  but  depends  on  the  course  of  the  computation. 
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In  particular,  the  adequacy  of  the  representation  at  each  step 
can  be  checked  by  using  one  more  term  and  verifying  that  Its 
effect  Is  small. 

In  order  to  apply  the  algorithm  just  described,  one  needs  an 

0  n 

Initial  function  f  ,  as   well  as  boundary  conditions  on  f  .   Care 

must  be  exercised  when  the  boundary  conditions  are  Imposed: 

f(u,u  )  at  a  boundary  may  be  Imposed  only  for  values  of  u,  u 

such  that  the  vector  (u,u  )  points  from  the  boundary  Into  the  gas. 

The  distribution  of  the  velocities   of  the  particles  coming  from  the 

fluid  and  hitting  the  boundary  depends  on  the   flow  and  cannot  be 

Imposed  arbitrarily  (see[7]).   If  this  obvious  condition  Is  not 

respected,  numerical  Instability   will  result. 

Application  to  a  shock  problem.   Consider  a  gas  of  elastic  spheres 

flowing  In  -°°  <  X  <  +«>,  with 


(11)       f(-co,u,u^)  =  p^7T-^U-2exp(-((u-v^)^+u2)/u^) 


(12)       f(+«',u,u^)  =  p^-n-^\J^^exp(-{{u-v^)^+u^)/\J^^) 


clearly  u(-«>)  -  v  ,  u(+«>)  =  v  . 


The  Mach  number  M  Is  defined  by 


(13)  M  =  /E75  v^/U^  . 


11 


If  M  >  0  a  shock  will  develop.   There  may   be  a  steady  shock  if 
the  following  conservation  laws  are  satisfied 

(I'i)  PjV^  =  P2V2 


(15)  Pi(v^+  |u2)  =  P2(v2+  |u^) 


(16)  P  v-^(v^+  |u2)  =  P^v./vl+   |u^)  , 


where  it  is  assumed  that  the  ratio  of  specific  heats  is 


Y  =  5/3 


From  equation  (14),  (15),  (I6)  we  may  deduce 


(17)  (U^/U^)^  =  (m2+3)(5M^-1)/16m2 


(18)  (v^/v^)   =(M^+3)/4M^ 


An  important  parameter  in  the  shock  problem  is  the  shock  thickness, 
conventionally  (and  awkwardly)  defined  by 


V  -v^ 
(19)  X  =  -^   ^ 


I  du  I 
max  :j— 
X  '  dx' 


12 


We  pick  p-^  =  1,  v-^  =  1.   Given  M,  equation  (13),  (17),  (l8) 
yield  U.^  ,  pp,Vp,  U  .   We  call  the  left  hand  end  of  the  shock 
the  upstream  side;  we  thus  chose  the  units  so  that  the  mean  free 
path  upstream  is  one.   In  those  units  X  is  the  ratio  of  the  shock 
strength  to  the  upstream  mean  free  path;  several  authors  have  studied 
X~   as  a  function  of  the  Mach  number  M. 

For  practical  reasons  we  replace  the  region  -°°  £  x  <_  +°° 
by  the  region  -a  <_  x  <_  a,  where  a  is  chosen  large  enough  so  that 
any  further  Increase  in  a  will  have  no  noticable  effect  on  the 
shock.   At  X  =  a  we  impose  the  boundary  condition 


f(a,u,u^)  =  P2TT  -'-U2^exp(-((u-V2)^  +  \xl,) /^\)    for  u  +  v^  <  0 


and  at  x  =  -a  we  impose  the  condition 


f(-a,u,u^)  =  TT  ^U-^^exp(-((u-l)^+u^)/U?)  for  u  +  1  >  0  . 


We  divide  [-a,+a]  into  K  -  1  segments,  with  a  spatial  increment 


Ax  =  2a/K  . 


Our  aim  it  to  obtain  the  steady  shock  profile  as  the  limit, 
when  the  time  tends  to  infinity,  of  an  unsteady  flow  starting  from 
an  initial  function  f  =  f(x,u,u  0).   This  initial  function  should 
be  chosen  so  that  the  steady  limit  is  achieved  as  fast  as  possible. 
We  first  tried  initial  function  f   resulting  from  an  approximate 


13 


solution  of  the  Boltzmann  equation.  In  particular  we  tried  the 
solution  of  the  Mott-Smlth  u   theory  [12].   This  turned  out  to  be 
a  very  poor  choice.   It  Is  clear  that  the  convergence  to  the  steady 
limit  Is  Inherently  slow-lf  we  use  K  points  across  the  shock,  and 
if  the  stability  condition  (9)  Is  respected.  It  takes  at  least  K 

steps  for  the  fastest  particles  to  cross  the  shock.   If  f   Is  the 

8  f 
Mott-Smlth  solution,  the  Initial  values  assumed  by  jrr   are  very 

small,  and  the  relaxation  to  equilibrium  takes  an  extremely  long 

time,  (showing,  by  the  way,  that  the  Mott-Smlth  solution  Is  not 

a  very  good  approximation  to  the  real  f).   In  addition,  some  odd  effects 

appear:   at  low  M  the  Mott-Smlth  theory  overestimates  the  shock 
width,  yet  with  Mott-Smlth  Initial  data  the  shock  at  first  appears 

to  widen;  this  effect  can  also  be  observed  in  the  work  of  Havlland 

[9]. 

After  considerable  experimentation,  it  was  found  that  an 
appropriate  f   is  the  one  which  corresponds  to  a  shock  of  zero  width 

[f(-«>,u,u  )   for  X  5  0 
(20)  f(x,u,u^,0)  =^ 

|f(+o°,u,u  )   for  X  >  0  . 


The  initial  f  given  by  (20)  is  particularly  appropriate  when  one 
tries  to  determine  the  shock  width  X  as  defined  by  (19).  X  is  a 
local  property  of  the  shock  center,  and  with  the  data  (20)  x 

9  f 
approaches  equilibrium  values  long  before  -^^r  becomes  close  to 

o  t 

zero. 

It  is  worth  noting  that  from  the  numerical  point  of  view  the 
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determination  of  the  shock  width  X  is  a  comparatively  difficult 
undertaking,  since  it  requires  high  accuracy  in  the  region  of 
fastest  variation  of  f.   In  a  variety  of  other  problems,  e.g. 
problems  involving  the  interaction  of  a  shock  with  a  boundary, 
the  choice  of  initial  data  is  less  critical  and  the  computation 
is  easier  to  carry  out. 

We  now  apply  the  method  outlined  earlier  to  the  study  of 

the  shock  wave.   There  is  a  considerable  number  of  numerical 

n  n 

parameters  to  be  chosen:   the  centers  v  and  scales  u   of  the 

X  X 

expansion  (3),  as  well  as  the  number  (L-, +1 )  (Lp+1)  on  nonzero 
terms;   the  size  2a  of  the  region  of  integration,  the  spatial 
increment  Ax,  the  time  step  At,  the  number  of  quadrature  points 
N^NpN-^Nii  in  each  evaluation  of  the  collision  integral  and  the 

number  N   of  points  at  which  f"^   is  evaluated  given  f  . 

n      n 
We  choose  v   and  u   as  follows: 

X        X 


n+1    — n.  s 

v    =  u  (x) 

X        ^  ' 


(21) 

u"-'^  =  y2RT^(x) 


X 


i.e.  we  expand  at  each  step  around  the  mean  velocity  at  the  pre- 
ceding step  and  using  a  scale  determined  by  the  temperature  at 
the  preceding  step.   u  ,  T^^,  are  given  by  (4)  and  (5).   This  choice 
is  not  the  only  reasonable  one,  and  will  be  further  discussed  below, 

The  width  2a  of  the  region  of  Integration  was  chosen  by  trial 
and  error,  generally  around  25  mean  free  paths.  Ax  is  chosen  small 
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enough  so  that  any  further  decrease  in  Ax  will  not  affect  the 

outcome  of  the  calculation.   We  proceed  as  follows:   We  evaluate 

4^  which  enters  the  definition  (19)  of  X  using  both  the  formula 
dx 

,   ,  du  ^   k+l~  k-1 

^"^"^^  dx      2Ax 

and  _ 

.       .  du  ^   k+1   k 

^^^^  dx  -     ^^ 


which  are  of  different  orders  in  Ax;  when  they  are  in  substantial 
agreement  Ax  can  be  considered  sufficiently  small.   It  was 
found  that  Ax  of  order  1  (i.e.  one  mean  free  path)  is  generally 
adequate;  under  these  circumstances,  X  evaluated  with  the  use  of 
(23)  is  a  more  reliable  estimate  of  the  true  X,  since  X  is  a  local 
property  of  the  shock  center  and  an  estimate  using  (23)  depends 
on  the  values  of  f  in  a  smaller  neighborhood. 

The  stability  condition  (9)  gives  a  good  estimate  of  the 
appropriate  value  of  At.   We  usually  choose  At  to  be  0.8  times 
the  maximum  value  allowed  by  (9).   Higher  values  of  At  may  give 
rise  to  instability  in  the  presence  of  temperature  overshoots 
while  lower  values  lengthen  the  computation  without  increasing  its 
accuracy.   At  this  point  we  have  to  introduce  an  additional 
numerical  parameter.   It  is  readily  seen  that  the  stability  of  the 
scheme  f'^'*"-'-  =  Af^  would  imply  the  stability  of  the  complete  scheme 
(8)  if  only  the  integrand  on  the  right  hand  side  had  compact  support 
This  last  condition  is  not  satisfied,  but  f  does  decrease  rapidly 
with  increasing  |u|,  |u  |,  so  that  one  might  assume  that  condition 
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(9)  is  sufficient  for  stability.   Numerical  experimentation  shows 

this  to  be  the  case  whenever  L-,_<  3  and  L„  <_  3.   However,  when 

L^  or  !„  is  larger,  the  range  of  u,  u  over  which  f  is  not 

negligible  increases,  and  it  is  necessary  to  truncate  the  support 

of  f.   This  can  be  done  by  setting  f  =  0  whenever 

|u|  >  V  +  ^Xu  , |u  I  >  A^u  ,  where  u   is  the  scale  of  the  expansion, 

X  X     r  X  X 

V   its  center, C  is  the  largest  root  of  Hj^(u)  =  0  and  A  Is  a  constant 
larger  than  1.   When  X  >  1  such  truncation  leads  to  no  decrease  in 
accuracy,  since  the  expansion  in  Hermite  polynomials  is  not  uniformly 
valid  in  u,  u  .   We  generally  chose  X  -  1.1. 

We  generally  took   L-,  =  Lp,  equal  to  an  integer  L.   Clearly  we 
must  have  L  <  N;  on  the  other  hand  if  N  were  much  larger  than  L, 
information  would  be  generated  and  immediately  discarded;  so  we 
generally  chose  N  =  I,  +  1  (L  even)  and  N  =  L  +  2  (L  odd).   The 
difference  between  the  odd  and  even  cases  is  due  to  programming 
consideration  and  is  of  no  particular  significance. 

This  leaves  open  the  choice  of  L,  the  number  of  Hermite  poly- 
nomials in  each  of  the  variables  u,  u  .   It  would  be  natural  to 

'   r 

choose  L  so  large  that  a. .  S  0  for  either  i  or  j  larger  than  L. 
It  turns  out  however  that  a. .  decays  much  more  slowly  with  i  and  j 
than  expected,  but  that  the  presence  of  the  higher  terms  in  the 
expansion  affects  but  little  the  computed  values  of  X  and  of  the 
density,  mean  velocity  and  temperature.   For  example,  at  Mach 
number  M  =  1.6,  a^^  near  the  center  of  the  shock  tends  to  the  steady 
value  a„2i  -  --^j  ^^^  within  computational  error  there  is  no  difference 
between  the  value  of  X  computed  with  L  =  4  and  the  value  computed 
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with  L  =  3 J  i.e.  neglecting  S-r^u'      ^^   does  appear  therefore  that 
the  lower  moments  of  f  are  almost  independent  of  the  higher 
moments,  a  result  both  surprising  and  natural.   It  also  appears 
that  the  assumption  underlying  Grad's  thirteen  moment  approximation 
[5],  namely  that  the  coefficients  of  the  Hermite  polynomials  of 
degree  greater  than  3  aa?e  small,  is  not  correct  in  itself  but 
could  lead  to  correct  answers.  We  made  runs  with  both   L  =  3  and 
L  =  ^.   It  must  be  added  that  although  the  values  of  X  do  not  seem 
to  depend  on  L  provided  L  _>  3j  when  M  <  2,  the  initial  rate  of 
change  in  X  does  depend  on  L.   This  is  probably  of  no 

physical  significance,  since  the  initial  data  are  wholly  unrealistic. 
The  relationship  between  our  method  and  Grad's  will  be  the  object 
of  further  investigation  elsewhere.   It  should  be  noted  that   when 
L  =  3  our  f  is  represented  by  8  coefficients  a..,  taking  into 
account  the  fact  that  by  symmetry  a. .  =  0  for  odd  j;  when  L  =  4 

"'-  (J 

our  f  is  represented  by  15  functions.   This  compares  with  5  functions 
for  the  one-dimensional  case  of  Grad's  expansion. 

N,  ,  Np,  N-,,  Nk  are  also  chosen  by  trial  and  error.   We  must 
have  N-,  >  L/2,  Nu    >   L/2,  so  that  the  highest  moments  of  f  used  enter 
the  collision  integral.   It  was  generally  found  that  with  L  =  3  or 
L  =  4,  the  choice  N,  =  Np  =  N_  =  N^  =  3,  i.e.  8I  Integration  points 
for  every  evaluation  of  the  collision  integral,  was  quite  adequate. 
The  fact  that  such  low  values  are  adequate  is  testimony  both  to  the 
power  of  Gaussian  quadrature  and  to  the  aptness  of  the  representation 
(3). 

The  existence  of  conservation  laws  affords  a  natural  check  on 


18 


accuracy,  since  no  exact  conservation  is  built  into  our  scheme. 
With  the  initial  data  (20),  and  with  a  large  enough,  the  mass, 
momentum  and  energy  in  the  shock  region  are  constant.   The  magni- 
tude of  the  numerically  induced  variations  in,  say,  the  mass 
provides  a  reasonable  indication  of  the  accuracy  of  the  computation. 

In  tables  I  and  II  we  display  the  relaxation  from  the  initial 
data  (20).   In  table  I  the  mean  velocity  is  tabulated  as  a  function 
of  X  for  low  values  of  t/At  and  at  Mach  number  2;  this  should  give 
a  qualitative  picture  of  the  behavior  of  the  numerical  process.   In 
table  II  the  instantaneous  value  X   of  the  reciprocal  of  the  shock 
width,  the  maximum  of  [-r-pl  ,  the  location  of  that  maximum,  and  the 
computed  total  mass  Q  in  the  shock  region,  are  tabulated  as 
functions  of  t/At  for  M  =  2.  It  is  seen  that  ||^|  does  not  decay 
to  zero  fast,  if  at  all,  and  that  X~   oscillates.   In  each  run  we 
therefore  estimated  the  range  of  values  assumed  by  X   ,  defined  as 
the  range  between  the  last  maximum  and  last  minimum  of  X   .   It  is 
not  clear  whether  the  oscillations  ever  die  out.   They  are  amplified 
if  the  width  of  the  region  of  integration  2a  is  chosen  too  small, 
but  they  can  no  longer  be  decreased  by  a  further  increase  in  2a. 
The  location  of  the  maximum  of  |^|  recedes  In  time,  showing  that 
upstream  convergence  is  slower  than  downstream.   Similar 
observations  were  made  by  Haviland  [  9  ].   Q,  the  total  mass,  is 
evaluated  by 

K 


(24)  Q  =  ^P  (iAx)Ax  , 

1=0 
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Table   I 
u  as   a   function   of  x   and  t 


M  =   2,    At   =    .413,    Ax   -- 

--   1.5 

X 

t/At   =   1 

t/At    =    4 

t/At   =   8 

12.75 

1.000 

1.000 

1.000 

11.25 

1.000 

1.000 

1.000 

-9.75 

1.000 

1.000 

1.000 

-8.25 

1.000 

1.000 

1.000 

-6.75 

1.000 

1.000 

.999 

-5.25 

1.000 

.999 

.999 

-3.75 

1.000 

.999 

.998 

-2.25 

1.000 

.995 

.961 

-0.75 

.979 

.873 

.726 

0.75 

.461 

.499 

.535 

2.25 

.437 

.458 

.478 

3-75 

.437 

.444 

.460 

5.25 

.437 

.438 

.451 

6.75 

.437 

.437 

.444 

8.25 

.437 

.437 

.440 

9-75 

.437 

.437 

.438 

11.25 

.437 

.437 

.437 

12.75 

.437 

.437 

.437 
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Table  II 

Relaxation  to  a  steady  shock 

M  = 

=  2,  At  =  .413,  Ax  =  1.5 

location  of 

t/At 

x-1 

max  g^       max  ^^ 

Q 

1 

.61 

.056        +  .75 

44.31 

2 

.56 

.069        -  .75 

44.29 

3 

•  50 

.088        -  .75 

44.29 

4 

.44 

.097        -  .75 

44.31 

5 

.38 

.099        -  .75 

44.33 

6 

.32 

.096        -  .75 

44.35 

7 

.27 

.087        -  .75 

44.37 

8 

.27 

.072        -  .75 

44.40 

9 

.28 

.057        -  .75 

44.45 

10 

.28 

.043        -2.25 

44.48 

11 

.28 

.049        -2.25 

44.51 

12 

.26 

.055        -2.25 

44.54 

40 

.20 

.056        -6.75 

44.30 

41 

.19 

.060        -6.75 

44.27 

42 

.19 

.061        -6.75 

44.25 

43 

.21 

.060        -6.75 

44.22 

44 

.23 

.057        -6.75 

44.19 

45 

.24 

.052        -6.75 

44.17 

46 

.25 

.045        -6.75 

44.15 

47 

.24 

.045        -8.25 

44.13 

48 

.23 

.051        -8.25 

44.11 

49 

.21 

.056        -8.25 

44.10 
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Table  III 


Structure  of  a  shock 


M  =2,  t  =  9.5192 


X 

u 

P 

T 

H/p 

-12.75 

1.000 

1.000 

.300 

-1.398 

-11.25 

.999 

1.000 

.300 

-1.398 

-  9.75 

.999 

1.000 

.300 

-1.398 

-  8.25 

.999 

1.000 

.300 

-1.398 

-   6.75 

.998 

1.001 

.301 

-1.403 

-  5.25 

.977 

1.020 

.316 

-1.467 

-   3.75 

.831 

1.165 

.431 

-1.787 

-   2.25 

.670 

1.371 

.622 

-2.270 

-      .75 

.660 

1.392 

.680 

-2.479 

+     .75 

.616 

1.587 

.660 

-^299 

2.25 

.546 

1.937 

.628 

-2.011 

3.75 

.501 

2.213 

.609 

-1.802 

5.25 

.485 

2.342 

.602 

-1.712 

6.75 

.479 

2.381 

.601 

-1.692 

8.25 

.477 

2.379 

.605 

-1.708 

9.75 

.472 

2.367 

.612 

-1.743 

11.25 

.467 

2.349 

.619 

-1.788 

12.75 

.437 

2.285 

.623 

-1.883 
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Table  IV 


Coefficients  a.  . 


Mach  number  =1.6,  x=L5,    t=10.76 


1    =    0 

1  =   1 

1    =    2 

1  =   3 

1  =   4 

j    =    0 

1.447 

-.0002 

.48 

-.04 

-.38 

J    =    2 

.24 

-.19 

-.005 

.05 

.01 

J    =    ^ 

.11 

.08 

.004 

-.02 

-.006 

a. .  =  0  for  odd  J 
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Table  V 
Reciprocal  shock  width  X~   as  a  function  of  Mach  number  M, 


M  =  1.^4 


L  =  3        L  =  i| 
12  to  .13 


Gilbarg 

and  Zierinj 

Paolucci   Mott-Smith     et  al. 


136 


116 


181 


M  =  1.6      .22  to  .2H       .22  to  .2H 


164 


238 


222 


M  =  1.8      .18  to  .21 


.205 


.28^1 


M  =  2.0      .19  to  .25   .23  to  .29      .381 


235 


32^1 
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It  is  seen  that  Q  varies  little;  whatever  variations  there  are  can 
be  ascribed  to  the  inaccuracy  of  the  formula  (24). 

In  table  III  we  display  the  structure  of  a  typical  shock.  The 
mean  velocity  u,  density  p  ,  temperature  T,  and  Boltzmann  H  divided 
by  p  ,  are  given  as  functions  of  x,  for  M  =  2  and  t  =  9. 5192.  The 
familiar  features  of  the  shock  appear:  u  and  p  vary  in  a  monotone 
fashlon;T  exhibits  an  overshoot,  see  [16];  H/p,  which  is  determined 
up  to  an  additive  constant,  displays  a  dip.  H  is  evaluated  from  f 
using,  as  usual,  Gauss-Hermlte  quadrature. 

In  table  IV  we  present  the  coefficients  a. .  for  x  =  1.5, 
M  =  1.6,  t  =  10.76.   The  purpose  of  the  table  is  to  show  that  a^^ 
at  that  point  is  not  small. 

Some  of  the  more  interesting  results  are  grouped  on  table  V, 
where  the  ranges  of  oscillation  of  X~   for  Mach  numbers  1.4,  1.6, 
1.8,  and  2.0  are  given,  with  both  L  =  3  and  L  =  4,  and  compared  with 
the  values  of  X   computed  by  Gilbarg  and  Paoluccl  using  the  Navier- 
Stokes  equations,  and  by  Mott-Smith  and  Ziering  et  al  using  their 
respective  theories.   As  expected,  at  M  =  1.4  the  computed  X~   is 
very  close  to  the  Navler-Stokes  result.   At  M  =  1.6,  where  the 
result  is  seen  to  be  independent  of  L  _>  3,  the  shock  is  thinner 
than  the  Navler-Stokes  shock,  with  X~   close  to  the  value  given  by 
Ziering  et  al.   Although  our  method  is  clearly  Inspired  by  Grad ' s 
work,  and  although  some  of  Grad's  ideas  are  resoundingly  vindicated, 
the  numerical  results  do  not  agree  with  Grad's,  whose  shocks  are 
always  thicker  than  the  Navler-Stokes  shocks.   It  seems  that  five 
moments  are  just  one  or  two  short  of  giving  an  accurate  description 
of  the  shock. 
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Between  M  =  1.6  and  M  =  1.8  there  seems  to  be  a  change  of 
regime;  suggestively  this  occurs  In  the  region  where  Grad's 
approximation  breaks  down.   Above  M  =  1.8  the  results  seem  to  agree 
with  the  Mott-Smlth  predictions. 

Comparison  of  these  results  with  available  Monte-Carlo  results 
is  difficult,  since  the  Monte-Carlo  calculations  in  the  literature 
cover  time  spans  too  short  to  be  of  any  significance.   The  results 
contradict  the  conclusions  of  Bird  [1],  whose  shocks  are  always 
thicker  than  the  Navier-Stokes  shocks,  and  they  are  in  some 
qualitative  agreement  with  the  conclusion  of  Haviland  [9],  but  one 
may  wonder  whether  this  is  more  than  coincidence. 

Generalizations  and  comments.  It  is  quite  clear  that  the  procedure 
of  the  preceding  section  will  break  down,  for  a  fixed  number  of 
terms  in  the  Hermite  expansion,  whenever  the  Mach  number  is  large 
enough;  certainly  by  the  time  all  the  velocities  u  =  v  +  u  ^., 
i  =  0,...,N,^^  roots  of  Hj^(u)  =  0,  are  of  the  same  sign.   With  u  , 
v  given  by  (21)  and  N  =  5  this  breakdown  occurs  just  above  Mach 
number  2.   One  could  keep  Increasing  the  number  of  polynomials  as 
M  increases;  it  is  more  reasonable  to  systematize  the  Mott-Smlth 
and  Ziering  et  al  procedures  by  representing  f  as  a  sum  of  two 
series  of  the  form  (3),  with  scales  and  centers  determined  respectively 
be  the  conditions  upstream  and  downstream  from  the  shock. 

Other  changes  in  the  scaling  (21)  may  be  justified:   for  example. 
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It  is  probably  beneficial  to  introduce  two  distinct  scalings  for 
the  variables  u  and  u  . 

Another  modification  our  basic  method  was  explained  in  [3]: 
the  evaluation  of  the  collision  integral  may  be  performed  by  Monte- 
Carlo  quadrature,  with  the  possible  help  of  the  variance  reduction 
technique  introduced  in  [3]-   This  should  be  particularly  effective 
close  to  equilibrium  when  the  integrand  of  the  collision  term  is 
small,  provided  this  term  is  not  separated  into  gain  and  loss  terms, 
as  was  done  by  Nordsieck  [13]. 

The  methods  of  this  paper  are  readily  generalized  to  problems 
in  more  dimensions  and  with  other  types  of  Interparticle  force. 
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Appendix 

In  the  appendix  we  list  the  program  used  to  obtain  the  results 
given  above.    it  should  be  borne  in  mind  that  this  program  was 
written,  not  in  view  of  minimizing  computing  time  per  run,  but 
rather  so  that  experimentation  and  change  are  as  easy  as  possible. 
Obvious  ways  in  which  the  running  time  could  be  shortened  are: 
better  exploitation  of  the  symmetry  f(x,u,u  )  =  f(x,u,-u^);  use  of 
the  recursion  relation  between  Hermite  polynomials,  and  restriction 
of  the  calculation  to  regions  in  space  where  substantial  changes 
are  occuring.   As  it  stands,  the  prograa  performs  one  time  step 
in  approximately  a  minute  (with  K  =  l8,  L  =  ^ ,  N,  =  N„  =  N^  =  N^  =  3 
N  =  5)  on  the  CDC  6600  computer. 
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PROGRAM  COLLIDE  (OUTPUT) 

COMMON  CA.CB,CC.CD»CE,PI,TPI,TIKE,C1,C2»C3,C4,C5, 
1  C6,C7,C8,C9,Ci:,ROOT(10),WElT(10       ) 
3  .  POOT(10)'PFIT(10) 

DIMENSION  UU (5,5,40 )#VV(5,5, 40 ),D( 5,9,40) 
1  ,TT(40),UX(40).BOL(40),TRL(40) 
1  ,URL(40) 

IPRs-1 
C  IPR  DETERMINES  WHETHER  TC  PRINT  OR  HlT 

C         TT  IS  THE  SQUARE  OF  THE  THERMAL  VELOCITY 
C  TRLsRHO*TT=D(l,l, IC)*TT(IC) 

C         RATIO  OF  SPECIFIC  HEATS  QAMMa=5/3 
M  =  3 
M  =  5 
MM  =  3 
MMMa3 

CALL  HERI1Y 
1  (M,MMM) 
PIH3PI/2. 
COsO. 

PRINT  9n2n,C0,Cl,C2,C3,C4 
NC0UN=5 
IC0UN=1 
EPSsO. 
TIMEsO. 
NC  =  18 

NHALF5NC/2 

NHALFPsNHALF*! 
NCM=NC-1 
CNCbNC-1 
DX=1.5 
XXX3DX*CNC 
C         LA, LB  NUMBERS  OF  POLYNOMIALS 
LA  =  4 
LB  =  LA 
C         RR  UNIVERSAL  GAS  CONSTANT 
RRsO.5 
XMACH=2. 

XMMsXMACH*XMACH 
Ul  =  l. 
Rl  =  l. 
CClaUl*Ul*(6./5.)/XMM 

U2bU1*( (XMM*3. )/(4.*XMM) ) 
CC2aCCl*((XMM*3. )*(5.*XMM-1. )/C16.*XM^)) 
R2aRl*Ul/U? 
E1*CC1*R1 
E2=CC2*R2 

PRINT  9025 

PRINT  9017,  XMACH,U1,R1,CC1,EI,U2#R2,CC2,E2 
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c 
c 
c 


c 
c 
c 
c 
c 
c 
c 
c 


UXMs 

UMAX 

1*1.2 

1  *ux 


COFB 
DIST 
TIME 


8AM  =  R1*(L'1*U1*0.5*CC1) 
BAMMsR2*(U2*U2*0.5*CC2) 
PRINT  9C20,  BAM,BAMM 


UXM  IS  THE  LARGEST  MEAN  VELOCITY 
1.2 
=ABS(ROOT(M) ) 

h 

DT=DX/UMAX 
PRINT  9023,rT,DX 
XLAMsDT/DX 
=DT/(SQRT(2. )*PI ) 
ANCE  MEASURFD  IN  MEAN  FREE  PATHS, 

IN  COLLISION  TIMES,  VEL  IN  TH  VEL 


UNITS 


C 
C 


SET  ALL  TO  ZERO 
DO  120  1C=1,NC 
BOL(  IC)  =  0. 

DO    120    lA8l,5 
DO    120    IBai.S 
D(  lA,  IB,  lOsO. 

120  CONTINUE 

SETTING  INITIAL  DATA 
IC  =  1 
ICsNC 

DO  121  ICsl^NHALF 

TRL( IC)sEl 

UX( IC)=U1 
URL( IC)3UX( IC) 

TT( IC)=CC1 

D(l,l, IC)«R1 

121  CONTINUE 

DO    122    ICaNHALFPiNC 
UX(  IC)aU2 
TRL(  IC)=E2 
URL(  IC)=UX( IC) 
TT( IC)=CC2 
Dd/l.  IC)  =  R2 

122  CONTINUE 


DO    731     ICad,NC 
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CD=2.*RR*TT( IC) 
CE=PI*CD 
DO  500  l3l,H 
DO  500  J=1,M 
OX  =  ROOT(  I) 

QY=ROOT(j) 
GLOGsO. 

DO  501  IA=1,LA 
IAM=IA-1 
DO  501  IB=1,LB 

IBM=IR-1 
GL0GbQLOG*D( IA*I8, IC) 
1*H(IAM,QX) 
2  *H( IBM,QY) 
501  CONTINUE 

UU( I»J, IC)=QLOG*EXP(-QX*GX-QY«QY)/CB 
500  CONTINUE 

IF(IPR.LE.O)  GO  TO  1B06 
PRINT  9004 
PRINT  9022/ IC 
PRINT  9004 
DO  790  Jsl,M 
790  PRINT  9001, <UU(I,J,IC)/I»1,M) 
1806  CONTINUE 
731  CONTINUE 


C 

c 
c 
c 
c 
c 
c 


c 
c 
c 
c 


TIME  STEP 
NSTPai2 
NSTPslo 

NSTP=20 
NSTPs35 
NSTPS50 

DO  300  ISTP=1,NSTP 
PRINT  9007, ISTP 
TIME  =  TIME-^DT 

PRINT  9026,  TIME 

A  POINT  IN  SPACE 


DO  700  lCs2,NCM 
Sin=RR 

CD32.*RR*TT(IC> 


TIME  STEP    TIME  STEP 
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c 
c 
c 
c 


CC=SORT(CD) 
CE=PI*CD 

CA  =  CC 
CB  =  CC*Cil 
CCT=CC 
RSCL=SQRT(TPLf IC )/(D{l»l#IC)*TT(lC))) 

CONVECTIO:-.  TERMS 


DO 


M 


)    IP: 
IP  =  1 
IP: 


IC-1 


NC 


730  1=1, 
IP=IC*1 
VELC=ROOT( I )*RSCL 
i  ♦C1*D{2,1, IC) 
1  /D(l,l, IC) 
VEL5VELC*CC*UX( IC) 
IFCVtL.GF.O 
1F( IP.LT.:  ) 
ir(  IP.G^  .NO 
CP=ABS(VFL»DT)/DX 
CENT=l.-rp 
RATP=TT( IC)/TT( IP) 
CEP=CE/RATP 
RATP=SQRT(RATP) 
VELPaVELC*(UX( ir)-UX(IP))/CC 

VELPcVELP*RATP 
DO  730  J«1,M 

VELRC=ROO^( J) 
1*RSCL 

VELRsVELRC 
VELRPsRUOK J)*RATP 
1*RSCL 

CAMEsO- 
DO    710     IAsl,LA 
lAMslA-l 
DO    710     IB=1,LB 
1         ,2 

IBM=IB-1 

*CENT*D( lA, IB, IC)*H( I  AM, VElC ) *H( IBM, VEL9C) 
*  EXP(-VELC*VELC-VEURC*VELRC)/CE 
*CP*D(IA,IB,IP^*H(IAM,VELP)*I-(IBH,VELRP) 
♦EXP(-VELP*VELP-VELRP*VELRP)/CEP 
CONTINUE 

VV(  I, J. IC)  =  CAME 
730  CONTINUE 

ir(  IPR.LE-O)   QO  TO  180C 
PRINT  9022, IC 


1 
2 

1 
2 


710 
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c 
c 
c 
c 
c 
c 
c 
c 
c 


c 
c 
c 
c 
c 


PRINT  9004 
PRINT  9020,C^,CB,CC/CD,CE 

PRINT  9004 
DO  734  J=1,M 
734  PRINT  9001, (VV(I,J,IC)»I=1#M) 
1800  CONTINUE 


VARIOUS  POINTS  IN  U-SP^CE 

DO  100  jQsl,K 
DO  3.00  10  =  1, M 
RSLT=0. 
UAPsROOT( IQ) 
1*RSCL 

1  +C1*D(2,1, IC) 
1  /D(l,l, IC) 

UBPsROOT( JO) 
1*RSCL 

CATAaO. 

DO  410   lAal,LA 
lAMsIA-1 
DO  410  lBsi,LB 

1    ,2 

IBMelB-l 

CATAsCATA*D(IA,IB, IC) 
1  *H(IAM,UAP)*H(IBH,UBP) 
410  CONTINUE 

CONSTANTS        PI*PI/4     (FROM  ANGLES) 
♦  CD(TWO  INTFGRATIONS)   CC  (FROM  VE  )  /  PUPI*CD*CD 
FROM  NORMALIZATION  OF  FF,  /2  FROM 
RESULTING  IN  DIVISION  BY 
8*CC 
QsEXP(-UAP*UAP-UBP*UBP) 
RUG=8.*CC 
0=Q/RUG 


CALL 


HERMY(MM 
DO  10 
DO  10  . 
DO  10 


MMM) 
JA=1,MM 
JB=1,MM 
Jf;al,MMM 


DO  10  JDsl,NMK 
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c 
c 
c 


AN 

B 

A 

C 

SA 

B8 

BB 

VE 

1  ^ 

1  * 

UB 
UC 
UB 

QA 
GA 
GA 

DO 


lA 


1  * 

2  * 
205  CO 

DO 
lA 


1  * 

1  * 

210  CO 


NE 


GLE 
BA 
B  =  S 
ZsP 
AZs 
Z  =  S 
C  =  B 
S  =  B 
=  (C 
(CO 
COL 
UA 
=  UB 
=  U6 
=  SQ 
G 
USB 

use 

USB 
D 

20 
D 

,2 
IBM 
M=l 
DAT 
H(  I 
H(  I 
NTI 

21 

M=I 

DO 

,2 

IBM 

PAT 

H(  I 

H(  I 

NTI 

Q 

D 

DO 


c 
c 

=p 

=  c 

IN 

IH 

CO 

IN 

B* 

B* 

OL 

LL 

LB 

=  U 

P* 

P  + 

RT 

AU 

rC 

=  C 

=  S 

AT 

5 

0 


COLL 

OLLA 

OLLB 

IH*P 

0S(A 

(ANG 

*PIH 

S(AZ 

(AZ) 

CAZ 

SAZ 

LA-U 

B*CA 

♦SAZ 

AP*V 

VE*B 

VE*B 

(UB* 

SA  =  C 

OLLB 

OLLB 

QRT( 

AsO. 

lAsl 

205 


iSior' 

=Rf 0T( JA) 

=RCO"( JB) 

IH*pnOT(  JO 

MGLE) 

LE) 

♦POOTC JD) 

) 


AP) 

Z-U 

♦BB 

£♦8 

BC 

BS 

UB  + 

OLL 

-VE 

-VF 

GAU 


*RA 

BP)*BBC 

S 

A 


ur*uC) 

A-VE*BA 

*BBC 

♦  BBS 

SB*GAUSB*GAUSC«GALSC) 


,LA 
IB=1.LB 


=  IB 

A-1 
A  =  D 
AM, 
BM, 
NUE 
PAT 
0  I 
A-1 
210  IB=1,LB 


ATA*D(  lA, IB, IC) 

UA) 

UB) 

AsO. 
A=1»LA 


=  18-1 

A  =  PATA-^D^  lA,  IB/  IC) 
AM,GAUSA) 
BM.GAUSB> 
NUE 

A  T  A  =  0  . 

0  411  1A=1-LA 
IAMaIA-1 
411  IB=1,LR 
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1    .2 

IBMsIB-1 
QATAsQATA*D( lA, IB# lO 
1  *H(  IAM,COLLA)*F  (IBM/COLUB) 
411  CONTINUE 

AUX5(DATA*PATA-(.ATA*CATA)*ABS(VE) 
1  *SIN(ANGLE) 

RSLTrRSI.T+AUX*WEIT(JA)«VvEIT(JB)*PElT{jC) 

1  ♦PEIT(JD) 
10  CONTINUE 

CALL  HERMY(M,C) 
30  CONTINUE 

IFCRSLT.LE.P. )  RSLTsQ. 
UUdO.jQ,  IC)sRSLT*COFB 
1  *Q 
1  ♦VVC IQ» JQ, IC) 

100  CONTINUE 

ir( IPR.LE.O)    QO  TO  1801 

PRINT  9016 
DO  302  J=1#M 
PRINT  9001, (UU(I#J#IC), 1=1, M) 

302  CONTINUE 

1801  CONTINUE 

700  CONTINUE 

c 
c 
c 
c 
c 
c 
c 
c 

C         RESCALING 

DO  1707  ICb2,NCH 

TT( IC)=TRL( IC)/D(1,1# IC) 

UX(IC)=URL( IC) 
1707  CONTINUE 
C 

c 

C  ASSUMES  RR*?.=1. 

IC0UN=IC0UN-^1 
IF( ICOUN.QT.NCOUN)     IC0UN=1 

c 
c 
c 
c 

DO  1700  ICsl.NC 
CD=2.*PR*TT(IC) 
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c 
c 
c 


c 
c 
c 


CC=SORT(CD) 
CA  =  CC 

cb  =  cc*c:d 


ANALYSIS  or  DATA 
DO  177  IA-1,LA 

DO  177  IB=1,LH 

D( lA, IB, IC)=0. 
177  CONTINUE 

DO  102    J=1.M 
DO  102  I=1»M 

WElGHsWni^( I  )*WEIT(J) 
1  *CA 

UUU=UU( I . J, IC) 
ZA  =  ROOT(  I  ) 

ZB=ROOT( J) 
AA  =  EXP{ZA«-ZA*ZR*ZB)*UUU*WEiGH*CA 

DO  102  IA=1,LA 
DO  102  lBsl,LB 
IAM=IA-1 
IBM=IB-1 

D(  lA, IB» IC)  =  D( lA, IB, IC) 
1  +AA*H( IAM,ZA)*H( IBM,ZB) 
102  CONTINUE 


ENTROPY 

RH0=D(1.1, IC) 
B0L( IC)=n. 
DO  750  Isl,M 
DO  750  J*1,M 
DO  750  Kxl,M 
ZA  =  ROOT(  I) 
ZB=R0OT(J) 
ZC=ROOT<K) 
ZR=SQRT(ZB*ZB*ZC*ZC) 

VAL=0. 

DO  751  lAsl,LA 
IAM=IA-1 
DO  751  IBsl,LR 
1  ,2 
IBM=IB-1 
VAL=VAL*D( lA, TB, IC)*H( IAM,ZA)*H( IBM,ZR) 

751  CONTINUE 

VALE*EXP(-ZA*ZA-ZR*ZR)/CE 

VALEsVAL*VALE 
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c 
c 

c 


IF(VALE.LE.l.E-6)  GO  TO  750 

VAL=VAL*ALOG(VALe) 

BOL( IC)=BOL( IC)+VAL*WEIT( I>*WEIT(J) 
1*WEIT(K) 
750  CONTINUE 

COTsRHO*PI 
1  /CC 

BOL( IC)rBOL( IC)/COT 


TEMPERATURE 
TREAL  =  CD*{ (3./2. )*D(1,1, IC)*SCRT(0.5)*D(3.1,  IC) 
1  ♦  SQRT(2.  )*D(1,3,  lO) 
TREAL=TREAL/3. 

TREAL=TREAL/RR 
UPDATING 

URL(IC)=UX(IC)fri*cC*D(2»l# lu) 
1  /D(l,l, IC) 

TRL(IC)=TrEAL 
IF(  ICOUN.ME.NCOUN)   GO  TO  764 
9022,  IC 


PRINT 
PRINT  9004 
PRINT  9020» 
PRINT  9004 

PRINT  9004 


CA,Ca,CC#CD,CE 


IB8a,LB 


PRINT  9015 
DO  232 

PRINT  9020, (D(IA, IB, IC),lA=l,LA) 
232  CONTINUE 

PRINT  9004 


764  CONTINUE 
1700  CONTINUE 


C 

c 
c 

c 
c 
c 


PRINT  90C4 


CHECK  ON  CONSERVATION 
PRINT  9010 
SAsO. 
SB  =  0. 
SC  =  0. 
ABRAsl. 

DO  221  IC=1,NC 
QZ=UX( IC) 
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*CUF) 


221 


C 
C 
C 
C 


Q7.  =  URL(  IC) 

QUQU=D(1.1, IC) 
QUP  =  0.5*TRL(  lO/QUQU 
QA=QU0U*Q7 
QB=QU0U*(Q7*QZ*QUP) 

QC  =  QU.':u*n7*{QZ*QZ*5. 
SAsSA*QA 
SB=SB+Q8 
SC  =  SC-^-Da»l'  K) 

PRINT  9018,  QA,QB,QC 

CONTINUE 
SA=SA*DX 
SB=SB*DX 
SC=SC*DX 
PRINT  9004 
PRINT  9004 
PRINT  90lfl,SA,SP,SC 
PRINT  9004 
PRINT  9004 
470  CONTINUE 


SUMMINU  UP 

PRINT  9004 

PRINT  90 C 8 

PRINT  9025 

DO  222  IC=1,NC 

STR=0 . 

HI=TRL( IC)/D( j ,1, IC) 

PRINT  90  2  0,URl.  (  IC).n(l»l#  IC)  ,FI»BOL(  IC) 
CONTINUE 


222 


C 
C 
C 
C 
C 


760 


761 


EXIT  AND  REFINEMENT 

ERR  =  0. 

DO  760  IC=1,NC 

CRT=ABS(URL( in)-UX(IC) ) 

ir(ERR.LE.CRT)   ERR«CRT 

ir(ERR.LE.CRT)   INQsiC 

CONTINUE 

ERR=ERR/DT 
PRINT  9002»F-RR 
PRINT  9006, iNG 
IF(ERR-EPS)   761,762.762 
CONTINUE 
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CALL  EXIT 
762  CONTINUE 
C 

c 

C         SHOCK  THICKMESS 

TDX=2.*DX 

SLOPE=0. 

SLOPC=0. 

DO  1900  ICs?l»NCH 

SLXC=ABS(URL( IC-1)-URL( IC-D )/TCX 

IF(SLOPC.LE.SIXC)   SUOPCsSLXC 

IOT-IC-1 

SLX=ABS(URL( IC ) - URL( IC-1 ) )/DX 

IF(SLOPE.LE.SLX>       SLOPE=SLX 
IFCSLOPE.LF.SLX)    INQ=IC 
1900  CONTINUE 

SHW=(U1-U2)/SL0PE 

SHWC=(U1-U2)/SL0PC 

PRINT  90  28.  SHW,  SHWC 

Q=1./SHW 

QG=1./SHWC 

PRINT  9003, G,GQ 

PRINT  9006, ING 
300  CONTINUE 
C  END  OF  TIME  STEP 

C 

c 
c 
c 

C  COMPUTATION  OF  NORMAL  VELOCITY  COMPONENT 

PRINT  9024 
C         LOX  MUST  BE  ODD 
L0X=21 

L0XH=L0X/2 

CL0XH=L0XH 

QAMUTS2. 

DN=GAMUT/CLOXH 

DO    770     IRT3l,3 

IF(IRT.EQ.l)     IC=1 

IF(IRT.EQ.2)     IC=NHALF 

IF(IRT.EQ.3)     IC=NC 

PRINT    9004 

PRINT    9022, IC 

PRINT    9004 

CC=SQRT(TT( IC)) 

CE=PI*TT( IC) 

DO    770     U1,L0X 

IM=I-1 
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c 
c 
c 
c 
c 
c 


CIMsIM 

PLCs-GAMUTt-niM^-.N 

PLC=PLC-UX(  C) 

PLC=PLC/CC 

CHAT  =  0  . 

DO  901  J=1,M 

PLCN=ROOT( J) 

CHIEU=0. 

DO  900  IA=1.LA 

IAM=IA-1 

DO  900  IB=1.LR 
1  .2 

IBM=IB-1 

CHIEN=CHIEN 
1  *D( lA, IB, IC)*H( IAM,PLC)*H( iBf ,FLCN) 

900  CONTINUE 
CHIEN=CHIEN*EXP( -PLC*PLC) 
CHlEN  =  CHIEN*CC/i-E 
CHAT=nHAT*CHI£N*WEIT(J) 

1*PI 

901  CONTINUE 

PRINT  9020,  CHA' 
770  CONTINUE 


9001  FORMAT  (lX,6Fll.7) 

9002  rORMAT(4H  ERR.  F11.7) 

9003  F0RMAT(12H  REGIT RQCALS, 2F11 . 7 ) 

9004  FORMAT(IX) 

9005  F0RMAT(6H  EXArT,F11.7) 

9006  F0RMAT(1X,2 U5) 

9007  F0RMAT(/5H  STEP, 15/) 

9008  F0RMAT(14H  SHOCK  PROFILE) 

9009  F0RMAT(5H  TEMP,  F11.7) 

9010  F0RMAT(13H  JONSFRVAT I  ON ) 

9011  FORMAT  (5H  RSLT . Fll . 7 , 3HVAR , Fll 

9012  F0RMAT(6H  POI\T, 15) 

9013  F0RMAT(16H  IS  T^ E  SOL  GOOD,  F11.7) 

9015  FORMATCSH  C')EFF  D) 

9016  F0RMAT(3H  UU/) 

9017  FORMATOH  MACH  NO,  Fll .  7/lX  ,  4F  11 .  7/lX .  4F11 . 7/ ) 

9018  F0RMAT(ix,3(Fll.7,5X) ) 

9020  F0RMAT(1X,5^"11.7) 

9021  FORMATdX, 4'  11.7,213) 


8) 
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9022  F0RMAT(6H  SPACE. 13) 

9023  F0RMAT(6H  DT  DX,  2F11.7) 

9024  rORMAT(14H  N  VEL  PROFILE) 

9025  FORMAT(lX,lHU.10X,lHR,10X,lHT.10X,3hBOL/> 

9026  F0RMAT{5H  TIME,  Fll.7) 

9027  F0RMAT(14H  NAVICR  STOKES) 

9028  FORMAT  (/12H  SHOCK  W IDTH, 2F11 . 7/ ) 
CALL  EXIT 

END 


SUBROUTINE  HERMY 

1  (M,N) 

COMMON  CA,CB,CC,CD.CE,PI,TPI,TIhE,Cl,C2»C3.C4,C5, 
1  C6,C7,C8,C9,Cin,ROOT(10),WElT(lQ       ) 
3  .  POOT(10).PEIT(10) 
9001  F0RMAT(9H  MMMMMMMM»  13) 

PI=3. 14159  26535  89793  23846 
TPI=2.*PI 
C 

C  GAUSSIAN  QUADRATURE 

POOT(l)sO. 

c 

C  7  POINT  QUADRATURE 

IF(N.NE.7)  GO  TO  7 

POOT(2)=0. ^058451513 

POOT(4)»0. 7415311855 
POOT(6)=0. 9491079123 

PEIT(l)s0. 4179591836 
PEIT(2)sO. 3818300505 
PEIT(4)=0. 2797053914 
PEIT(6)=0. 1294849661 
7  CONTINUE 

C 

C         5  POINT  QUADRATURE 
IF(N.NE.5)  GO  TO  5 

POOT(2)=0. 5384693101 

POOT(4)=0. 9061798459 

PEIT(l)30.56888888a8e 88 8888 88888 98 
PEIT(2)aO. 4786286704 
PEIT(4)=0. 2369268850 
5  CONTINUE 
C 

C         3  POINT  QUADRATURE 
IF(N.NE.3)  GO  TO  33 

POOT(2)»0. 7745966692 
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c 
c 


c 
c 


PElT(l)=0.e  8888  888<B8888  88  86  88  88  8 
PEIT(2)=0 .55555555555555555 
33  CONTINUE 

P00T(3)- -P00T(2) 
P00T(5»=  P00T(4) 
P00T(7)--pooT(6) 
PEIT(3)=PFIT(2) 

PeiT(5)=PEIT(4) 

PEITf 7)aPEIT(6) 

GAUSS  HERMITE  QUADRATuFE 
ROOT(1)^0. 
7  POINT  QUADRATURE 
IF(M.NE.7)  GO  TO  77 

RnOT(2)=0  8162878828 
R00T(4)s.L. 673551628 
R00T(6)=2. 651961356 

WEIT  (1)=0  8102646175 
WEIT(2)=0. 4256072526 
WEIT{4)=0. 05451558281 
WEIT(6)=0 . 000971781245G 
77  CONTINUE 

5  POINT  QUADRATURE 
IF(M.NE.5)  GO  TO  55 

R00T(2)=0. 9585724646 
R00T(4)s2. 02018287 
WEIT(1)=0  .9453087204 

WElT(2)aO. 3936193231 
WEIT(4)=0. 01995324205 
55  CONTINUE 


lr(M.NE.3)   GO  TC  333 
R00T(2isl. 224744871 
WEIT(1)=1. 1816359 

WEIT(2}30. 2954089751 

333  CONTINUE 

R00T(3)=-R00T(2) 
R0OT(5)=-RO0T(4) 
R0OT{7)»-R0OT(6) 
WEIT(3)=WEIT(2) 

WEITt5)aHEIT(4) 

WEIT(7)-WEIT(6) 


C 
C 


NORMALIZATION  FOR 
SQP=1. 


HERMITE  POLYNCMIALS 
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CD.CE 

ME 

C4,C5 

C7,C8,C9,C10 


C1=SQRT(2.*SQP) 
C1=C1/^. 

C2=SQRT(8.*SQP) 
C3  =  SQR"^(48  *SQP) 
C4=SQRT(384.*S0P) 

C5=SQRT(3840. ) 
C6=C5*SQRT(12. ) 
C7sC('>*SQftT(14.  ) 
C8=C7*SQRT(16. ) 
RETURN 
END 
FUNCTION  H{ I»X) 
COMMON  CA,CB,CC 
COMMON  PT,TPI .1 
COMMON  C1,C2,C3 
COMMOM  C6. 
II=I*i 

GO  TO  a, 2, 3, 4. 5, 6, 7. 8, 9)   II 
H  =  l. 
RETURN 
CONTINUE 
H=X/C1 
RETURN 
CONTINUE 

H2=4.*X*X-2. 
H=H2/C2 
RETURN 
CONTINUE 

H3=(a.t«'X*X-12 
H=H3/C3 
RETURN 
CONTINUE 

Y  =  X*X 
Hs(16 

H=H/C4 
RETURN 
CONTINUE 

Y  =  X*X 
H=( (32.*Y-160 
H=H/C5 

GOTO  10 
CONTINUE 

Y  =  X*X 

l|=((64.*  -480.  )*Y*720.  )*Y-120 

H=H/C6 
GOTO  10 
CONTINUE 
CONTINUE 


♦X 


*Y-48. )*Y+12. 


)»Y+120. )*X 
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10 


H  =  0. 

CONTINUE 
RETURN 
END 


EOF 
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This  report  was  prepared  as  an  account  of 
Government  sponsored  work.   Neither  the 
United  States,  nor  the  Commission,  nor  any 
person  acting  on  behalf  of  the  Commission: 

A.  Makes  any  warranty  or  representation, 
express  or  implied,  with  respect  to  the 
accuracy,  completeness,  or  usefulness  of 
the  information  contained  in  this  report, 
or  that  the  use  of  any  Information, 
apparatus,  method,  or  process  disclosed 
in  this  report  may  not  infringe  privately 
owned  rights ;  or 

B.  Assumes  any  liabilities  with  respect  to 
the  use  of,  or  for  damages  resulting  from 
the  use  of  any  information,  apparatus, 
method,  or  process  disclosed  in  this 
report . 

As  used  in  the  above,  "person  acting  on  behalf 
of  the  Commission"  includes  any  employee  or 
contractor  of  the  Commission,  or  employee  of 
such  contractor,  to  the  extent  that  such  em- 
ployee or  contractor  of  the  Commission,  or 
employee  of  such  contractor  prepares,  dis- 
seminates, or  provides  access  to,  any  infor- 
mation pursuant  to  his  employment  or  contract 
with  the  Commission,  or  his  employment  with 
such  contractor. 
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